SUBROUTINE powell(p,xi,ftol,iter,fret)
  
  USE nrtype; USE nrutil, ONLY : assert_eq,nrerror
  USE nr, ONLY : linmin
  USE sfutils, ONLY : func
  IMPLICIT NONE
  REAL(SP), DIMENSION(:), INTENT(INOUT) :: p
  REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: xi
  INTEGER(I4B), INTENT(OUT) :: iter
  REAL(SP), INTENT(IN) :: ftol
  REAL(SP), INTENT(OUT) :: fret
   INTEGER(I4B), PARAMETER :: ITMAX=30
  REAL(SP), PARAMETER :: TINY=1.0e-25_sp
  INTEGER(I4B) :: i,ibig,n
  REAL(SP) :: del,fp,fptt,t
  REAL(SP), DIMENSION(size(p)) :: pt,ptt,xit

  n=assert_eq(size(p),size(xi,1),size(xi,2),'powell')
  fret=func(p)
  pt(:)=p(:)
  iter=0

  do
     iter=iter+1
     fp=fret
     ibig=0
     del=0.0
     do i=1,n
        xit(:)=xi(:,i)
        fptt=fret
        call linmin(p,xit,fret)
        if (fptt-fret > del) then
           del=fptt-fret
           ibig=i
        end if
     end do
     if (2.0_sp*(fp-fret) <= ftol*(abs(fp)+abs(fret))+TINY) RETURN
     if (iter == ITMAX) call &
          nrerror('powell exceeding maximum iterations')
     ptt(:)=2.0_sp*p(:)-pt(:)
     xit(:)=p(:)-pt(:)
     pt(:)=p(:)
     fptt=func(ptt)
     if (fptt >= fp) cycle
     t=2.0_sp*(fp-2.0_sp*fret+fptt)*(fp-fret-del)**2-del*(fp-fptt)**2
     if (t >= 0.0) cycle
     call linmin(p,xit,fret)
     xi(:,ibig)=xi(:,n)
     xi(:,n)=xit(:)
  end do

END SUBROUTINE powell
